Seismic inversion is a process of estimating from seismic data a model of subsurface formation properties, such as its reflectivity or acoustical impedance. Typically, a convolutional model is assumed in which the data are given by the seismic wavelet convolved with a reflectivity series. A problem for inversion is that the seismic wavelet is unknown, and in addition, both the data and wavelet are band limited. The present invention overcomes these problems by first measuring or estimating the signature of the source, i.e. the seismic wavelet before any attenuation has occurred, and then solving simultaneously for both reflection or rock properties and for amplitude attenuation and velocity dispersion, caused by propagation through the earth from the source to the reflectors and to the receiver. The present invention overcomes problems with wavelet uncertainty, time varying or changing seismic wavelets, the presence of noise in the inversion window, and lack of low frequencies. In addition, the present invention estimates the attenuation quality factor of the earth, which can provide further information about subsurface rock and fluid properties.
The Convolutional Model
Seismic inversion is a process of estimating a model of subsurface properties, such as the reflectivity or acoustical impedance, from seismic data. Inversion may be performed using raw data or processed seismic data, in which the processing may include migration or stacking. Typically, the model is assumed to be a convolutional model, which assumes a linear rock-physics relationship between stress and strain and Newton's laws of mechanics. The convolutional model is implied in most seismic processing and interpretation (See R. E. Sheriff, Encyclopedic Dictionary of Applied Geophysics, 4th edition, p. 67).
The convolutional model is expressed as a seismic trace p(t) being a convolution of an embedded (equivalent wavelet) w(t) with a reflectivity function r(t), plus noise n(t).p(t)=w(t)*r(t)+n(t),  (1)where the symbol * donates the convolution operator. The reflectivity r(t) typically represents the change in impedance at the edge of each bed or layer in the subsurface as a function of two-way travel time, but it may include different rock physics parameters. A reflectivity series can be estimated from well logs, and then Eq. 1 can be used to derive a synthetic seismic trace, called the seismic-well tie, using an estimate of the seismic wavelet. The simple inverse operation of dividing the data in frequency domain by the convolution wavelet is called deconvolution. Amplitude losses or attenuation is generally ignored in the model and in the inversion. Attenuation causes considerable loss of high frequencies, and thus a wavelet that changes with greater propagation distance and longer times. Consequently, the inversion is usually performed using small time windows, in which the seismic wavelet appears stationary or unchanging. The seismic data may also be preprocessed to partially compensate for attenuation effects, but since the attenuation is unknown, this compensation is not accurate.
It is well known that this practice of using Eq. 1 for inversion to obtain the reflectivity has a number of problems. The chief problem is that the wavelet, particularly its phase, is unknown. In addition, both the data traces p(t) and the wavelet w(t) are band limited, but the desired reflectivity is not. Because both reflectivity and the wavelet are unknown, assumptions are generally made about the reflectivity, the most common being that it is “white” and has no frequency dependence. These problems are compounded by the fact that what is desired is not the reflectivity, i.e. the change in seismic impedance at interfaces, but its integral or the formation seismic impedance. Seismic impedance can be estimated by dividing (or inverting) the seismic trace by the wavelet and integrating over time, but missing low-frequencies are especially problematic.
Problems with missing low frequencies are illustrated in FIGS. 1A-1F. As shown in FIG. 1A, a source signature is assumed with an amplitude spectrum 101 that is missing low and high frequencies. An earth model including attenuation is assumed and a synthetic seismic trace is computed to give 103 (FIG. 1C). Dividing by the wavelet, and integrating the seismic trace gives the integrated trace 104 (FIG. 1D), which is a poor estimate of the true earth impedance, which is shown as 105. Seismic sources typically lack low-frequencies down to 0.1 Hz, but if we assume such a source with a source signature given by 102 (FIG. 1B) and the corresponding broad-band seismic trace 106 (FIG. 1E), then the integrated trace 107 better matches the seismic impedance 105. It differs because the trace still includes attenuation effects. It should be noted that there is no commercial seismic source that can generate the broadband signature shown in 102. The desire for such a source is well known, and this area is an active area of research. Making an active source with such low-frequencies is very difficult. The present invention is a way to overcome the need for such sources.
In order to obtain a realistic impedance section from narrow-band data, Oldenburg, et al. (Geophysics 48, 1318-1337 (1983)) made a different assumption to solve the inversion problem. He proposed what is today called a sparse-spike inversion in which the reflectivity is assumed to be made up of a spikes or delta functions for a small number N layers.
                              p          ⁡                      (            t            )                          =                              w            ⁡                          (              t              )                                *                                    ∑                              k                =                1                            N                        ⁢                                                  ⁢                                          z                k                            ⁢                              δ                ⁡                                  (                                      t                    -                                          τ                      k                                                        )                                                                                        (        2        )            
In this equation and subsequent equations, noise is ignored. Sparse spike inversion is illustrated in FIG. 2. This inversion is nonlinear and must be solved iteratively or recursively using model optimization methods. In this case, the unknown model parameters (zk and tk, plotted in 203) as well as the wavelet w(t) (202) are determined that minimize the misfit between the observed data pobs(t) and the model synthetic pmod(t) (201) to some power L.∥pobs(t)−pmod(t)∥L=minimum  (3)
The value L=1, the absolute-value norm, is preferred over the least-square norm L=2. It is well known that using the L=1 norm (or “L1” norm) yields a sparse or small number of events, and inversion with the L1 norm is called sparse inversion. When the wavelet is known and is completely stationary, i.e. for model data, then the full impedance can be obtained from the inversion without the need for low frequencies as shown by Oldenburg. Effectively, the reflectivity is assumed to be independent of frequency and the wavelet is stationary and contains all the frequency variation. This assumption is not valid for field data.
Sparse-spike inversion is currently one of the most rigorous methods for post-migration inversion that results in blocky realistic impedance sections, and a version is available in a number of commercial inversion software packages. It is known to be able to partially extend the band to lower frequencies for field data, but most often it is run in a band-limited mode or with a low-frequency model supplied by the interpreter or from well logs (see for example, Volterrani, U.S. Pat. No. 7,719,923). Sparse-spike inversion is usually performed post-migration as a 1-D inversion on stacked data or as a multiple-trace AVO-type inversion on prestack data as described by Routh, et al. (U.S. Pat. No. 7,072,767).
The fact that the wavelet is unknown limits inversion methods using the convolutional model and sparse-spike inversion in particular. A number of methods attempt to solve for the wavelet along with the reflectivity as described by Routh, et al. (U.S. Pat. No. 7,072,767), but constraints and assumptions are required. Oldenburg et al. (Geophysics 46, 1528-1542 (1981)) and others show that the solution is non-unique. More commonly, the wavelet amplitude spectrum is estimated from the seismic data, and the well-seismic tie is used to estimate its phase. A small variation in wavelet amplitude at the low-frequency edge of the band has a large effect on the phase of the data, especially for minimum-phase data, and accurately estimating the phase of the wavelet is very difficult. Thus, well data is desired, but even with well data, uncertainty regarding the wavelet remains.
Problems with the inversion also arise from the stationary assumption for the wavelet. This assumption means that the inversion can only be performed in isolated small time windows. If there are two targets at different times, then the inversion is performed separately for the two targets, and the output merged. The presence of noise, such as multiples or surface waves, within the window can compromise the results, because these events have different raypaths than the primary reflections and thus different effective wavelets. A more fundamental problem is that a small time window corresponds to poor resolution in frequency. This lack of frequency resolution will limit the effectiveness of the spike assumption and the ability to capitalize on the expectation that reflection coefficients are frequency independent.
Source-Signature Deconvolution
In its most general form, the convolutional model relates the response at a receiver p(t) to the source signature s(t) convolved with the impulse response of the earth g(t), normally called the Green's function, plus noise n(t).p(t)=s(t)*g(t)+n(t)  (4)
Source signature deconvolution is the process of removing the source signature, which must be known. It cannot be completely removed, since it is bandlimited, but it can be replaced with a more desirable bandlimited impulse response. For example, the signature from marine airgun sources is complex and ringy, because of oscillations of the air bubble and ghost effects from the ocean surface. A deep-water signature is sometimes recorded using a deep tow hydrophone or the signature is determined by modeling the response of an array of guns, each with a measured signature. The signature is then use to deconvolve the data to remove the effects of the bubble and ghost to generate data as if it had been recorded with a source with a more desirable pulse. The deconvolution filter is applied to the complete trace. It may also be used to shape the data to a desired zero-phase pulse.
Use of Source-Signatures in Inversion
Ziolkowski (“Why don't we measure seismic signatures?” Geophysics 56, 190-201 (1991)) proposed measuring the source signature citing problems with the traditional use of the convolutional model (Eq. 1). He concludes that the convolutional model cannot fit the reality of a changing convolutional wavelet in which the signature changes as it propagates through an absorbing earth. He points out that many geophysicists are unaware of possible ways to determine source signatures and discussed possible methods. A specific method to use the source signature for inversion of reflectivity was not provided.
Equation 4 can be rewritten asp(t)=s(t)*(attenuation(t))*r(t)+n(t),  (5)where now the response of the earth (the Green's Function) is decomposed into the reflectivity function as in Eq. 1 and propagation effects due to both absorption and interbedded multiples. (See Osman and Robinson, introduction in Geophysics reprint series No. 18, “Seismic Source Signature Estimation and Measurement,” pp. 1-14 (1996).). It can be recognized that from Eq. 5 and Eq. 1 that the inversion wavelet is given byw(t)=s(t)*(attenuation(t)).  (6)
The effective attenuation in Eq. 6 is ray-path dependent and changes along the ray-path and with distance travelled. Thus, the wavelet changes with ray-path distance and with time, requiring use of small time windows in inversion. For example, the effective wavelet will be different for different offsets used for AVO inversion methods because of increased attenuation over the longer offsets compared to the nearer offsets. Lazaratos (U.S. Pat. No. 6,516,275) teaches a method to spectrally balance the near and far offset seismic data to account for simple attenuation effects and for the processing filter of NMO stretch. Lazaratos determines attenuation compensation, assuming a single effective Q value, which along with the NMO-stretch compensation makes the spectrum of the near and offset traces match. Once the data are spectrally compensated, then a single wavelet can be used to perform AVO inversion.
Equation 6 is not used to derive a wavelet for inversion for reflectivity. In addition to problems with measuring the source signature, attenuation is unknown and difficult to estimate from surface seismic data. Attenuation, or Q, is most often measured by differences in the amplitude spectra measured by geophones at two depths in the earth, for example by a Vertical Seismic Profile. One of the problems for surface seismic data is there is no reference provided by measurements at two different depths from which a change in attenuation is determined Another problem is amplitude spectra are also affected by the interference of multiple events in a time window. Q estimation from surface seismic data has been largely limited to characterizing anomalous high attenuation regions in shallow sections, for example using Q tomography to estimate attenuation from shallow gas pockets (Hu, US Patent Application 2011/0273961). Accurate measurements of attenuation for reservoir zones could be used to improve processing and attenuation compensation. In addition, measurements of attenuation could provide useful information about reservoir rock and fluid properties. For example, attenuation may be related to the permeability of reservoir formations.
Accordingly, there is a need for an improved seismic inversion method, one that avoids problems with wavelet uncertainty and changes in the wavelet with time. The improved seismic inversion method should also be able to: invert the full trace without relying on small time windows, be robust with respect to noise or other events recorded at the same time as the primary reflections, have less need for well information, and potentially recover some of the low-frequency information. In addition, a measure for attenuation would be beneficial both to improve processing, i.e. compensation for attenuation, or to provide further information about the subsurface formations.